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Abstract 

We introduce two novel quantum Monte Carlo methods and employ them to 
study the superfluid-insulator transition in a two-dimensional system of hard- 
core bosons. One of the methods is appropriate for zero temperature and is 
based upon Green's function Monte Carlo; the other is a finite-temperature 
world-line cluster algorithm. In each case we find that the dynamical exponent 

is consistent with the theoretical prediction of z = 2 by Fisher and co-workers. 
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Recently, the study of the superfluid-insulator phase transition in systems of disordered 
interacting bosons has attracted considerable attention. As discussed by Fisher and co- 
workers , this transition is associated with a zero-temperature transition from a superffuid 
to an insulating state as the strength of the disorder is increased. In addition to experimen- 
tal studies in systems such as 4 He adsorbed in porous media, random magnets, Josephson 
junction arrays, exciton lifetimes in quantum well structures, vortices in bulk type-II su- 
perconductors, a huge amount of numerical work has been performed on models of these 
systems ||. 

Computer simulations are a natural candidate for studying this transition. The results 
of the simulations to date, however, are somewhat surprising as they seem to indicate that 
slightly different models lead to different values of the relevant critical exponents. The most 
significant disagreement exists between the numerical work on the Vallain form of the disor- 
dered XY model which has produced a value of the dynamical exponent z consistent 
with the original prediction of z = 2 given by Fisher and coworkers, and work on the hard- 
core boson model ||, which obtained z = 0.5. This difference was suggested in the latter 
work to be the result of (unspecified) low-lying excitations that are not encompassed in the 
models studied by others. Recently, these claims have been questioned [@,(|||- Most of the 
numerical work on these dirty-boson systems has been based upon finite-temperature quan- 
tum Monte Carlo (MC) methods, although exact diagonalization and large-N methods 
H have also been employed. 

The aim of this Letter is to introduce novel zero- and finite-temperature quantum MC 
methods and to use them to study the superfluid-insulator transition in a hard-core dirty- 
boson system. As we discuss below, these two methods have important, distinct conceptual 
and practical advantages over the standard algorithms. Our zero-temperature algorithm is 
based on the Green's function Monte Carlo (GFMC) method || and allows one to compute 
ground-state properties directly and to perform scaling analysis in the spatial dimensions 
only. The finite-temperature algorithm employed is a new way to implement the world-line 
algorithm that can reduce the equilibration and auto-correlation times by several orders of 
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magnitude. 

Our Hamiltonian is the standard hard-core boson model with random disorder: 

H = -JY.ib\b 3 +b)b l ) + Y.V l n 2 l (1) 

{ij) » 

where J is the overlap integral between neighboring bosons, (ij) denotes a nearest neighbor 
pair of lattice sites, b\ and b^ are the boson creation and destruction operators at lattice site 
i, and Hi = b\bi is the boson number operator which is restricted to have the eigenvalues of 
or 1. Vi is the on-site disorder and is assumed to be randomly and uniformly distributed 
over the interval (—V,V). We consider iV& bosons on a two-dimensional square lattice of 
L x L (= N) sites, with periodic boundary conditions in both spatial dimensions. In this 
work we restrict ourselves to half-filling (p = N^/N = 1/2). 

Two- dimensions occupies a special place in the theory of this model because of the 
predicted existence of a finite universal value of the dc conductivity |l| . Scaling theory also 
predicts that near the critical value of the disorder V c the singular part of the total (normal 
plus superfluid) compressibility k is expected to behave as where 5 = V — V c . If 

z = d = 2, then near V c this quantity is finite and independent of system size. We observed 
this behavior in our zero-temperature simulations. 

To apply the GFMC method [IDf, we define a kernel K = c — H . In an occupation 



number basis R, all matrix elements K(R, R') are readily obtained, and the constant c is 
chosen such that these elements are non-negative. GFMC projects out the ground-state 
wavefunction ip by iterating the equation ^( t+1 ) = Kip® from any positive initial function 
^(o) = ipj*. We can rewrite the equation in a form more appropriate for the MC process: 

ft t+1 \R) ocY,K(R,R')4> (t) (R'), (2) 

R' 

where tp {t) (R) = ip T (R)<ip® (R) and K(R,R') = ij) T (R)K(R,R^ l {R'). The trial wave 
function ?/y for the ground state serves as an importance function. In the MC process, v/j® 
at each t is sampled by a finite ensemble of configurations (walkers) each carrying a weight. 
A walker at R' at time t advances to time t + 1 by sampling R from K(R, R') / J2r K(R, R') , 



and determining a new weight from the product of the old weight and J2r K(R, R') . Splitting 
and combining procedures are applied at times to control the weights. A good i/jt can greatly 
reduce the fluctuations in the weight factors. 

Once the iteration has reached the asymptotic limit, ground-state properties can be 
calculated from the distribution of configurations. The projection of the ground state here 
is free of the Trotter approximation or other systematic errors. The ground-state energy 
is given by E = (HipT /Vt) • The compressibility k = (N^d 2 E / dN 2 )' 1 can be computed 
directly from E(Nb) by finite differences. 

In GFMC, calculations of quantities other than the energy can usually be performed 
only approximately because the result of the simulation is a distribution of configurations 
proportional to ipTipo- However, the superfluid density is a rather special property, being 
directly associated with the response of the system to moving boundaries: it is simply 
proportional to the long imaginary-time diffusion distance of the world-lines. By analogy 



with the arguments in jTT], a diffusion displacement D can be defined: 
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I[£(r? + '>-r?>)], (3) 



where boundary crossings in the x- and y-directions are included. It can then be shown that 
the superfluid density p s is given by 

c-E <D 2 (r)> 
Ps = —ry- lim 4 

With GFMC, we are essentially examining the diffusion over a finite portion of the (infinite) 
zero-temperature paths. In the full path integral, the density of configurations would be 
proportional to the ground-state wave function ip at each end of this segment. We replace 
this condition at one end by an approximate one in which the density is proportional to 
ipT instead of tpo. However, we emphasize that in the limit (r — > oo) that determines 
the superfluid density, these edge conditions have no effect, and the method is exact. We 
expect this method of computing the superfluid density to be quite valuable in many other 
applications, including the continuum. 



To extract the superfTuid density, it is necessary to extrapolate diffusion distance to 
infinite r. Due to the local motion of the paths, even non-superfluid states will propagate a 
finite diffusion distance at large r. In order to project out this localized effect, we use the 
following to model the asymptotic diffusion 

(D 2 (r)) /r = a' + (b'/r)[l - exp(-cV)], (5) 

where the constant a' determines the superfTuid density and b' and d are associated with 
"local" diffusion. The effects of higher ( > 2) orders in 1/t are small and different functional 
forms with the same asymptotic expansion produce identical answers. 
Our ipx is the product of a one-body term and a Jastrow factor: 

MR) = $1 exp[- £ f(r mn )}, (6) 

m,n=l 

where $i = exp(— aJ2i n iVi ~ bJ2<ij> n iVj)- The Jastrow factor contains long-range two- 
body correlations of the form 

f a, if r mn = 1 

f(r mn ) = , (7) 

[P/ r mn> otherwise 

where r mn denotes the distance between bosons m and n. The variational parameters, a, 
b, a, j3, and T], are optimized by minimizing the variance of the local energy. Due to the 
discrete nature of the problem, importance sampling with the long-range correlation in ip? 
can be implemented in a very efficient way. We have verified that our results, in particular 
the superfTuid density, remain unchanged except in statistical errors when these parameters 
deviate from their optimal values. 

The results of our zero-temperature simulations are summarized in Figs, p] and 0. We 
calculated the energy, compressibility, and superfTuid density for system sizes of 4 x 4 up to 
12 x f 2. The disorder magnitude V is the controlling parameter of the phase transition, and 
for each V, about f 00 realizations of the disorder are studied. We verified that our results 
for the 4x4 system are in excellent agreement with the exact diagonalization studies [0. 
Figure [I| presents the superfluid density as a function of disorder V for various system sizes. 
The transition is evident by the changes in p s with system size for various V. 
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To determine the precise scaling behavior, we used the following scaling relation for the 
superfluid density: 

p s (6,L) = L- z Y(L^5) (8) 

where Y is the scaling function. From this form, we see that when scaled with the correct 
value of z, curves for various L will have a point of common intersection that marks the 
critical point V c (5 = 0). Results for z = 0.5, 2.0, and 3.0 are shown in Fig. |2|. Clearly, z = 2 
is much favored over the other two values. A statistical analysis yields V c = 9.9 ± 0.4, and 
z = 2.0 ±0.4. As noted previously, if the z = 2 scaling behavior is correct, the compressibility 
k should be finite and continuous in the region of the transition. We observed no statistically 
significant change in k in this region (see Fig. [I] inset), and found k ~ 0.10(1). We also 
determined the critical exponent v to be 0.9 ± 0.1. This result is consistent with the lower 
bound of unity predicted by theory and with the value obtained for the quantum rotor 
model. It is not consistent with the estimate of 2.2 ± 0.2 found in |J. 

To obtain evidence complementary to the zero-temperature calculation, we also per- 
formed finite-temperature MC simulations and employed a finite-size scaling analysis with 
two lengths, the lattice size L and the inverse-temperature (3. 



To produce data, we extended a recently developed cluster algorithm [TJj which is a 
form of world-line MC found to have smaller equilibration and autocorrelation times than 
the conventional implementation. The algorithm is based on two simple observations: in 
the absence of disorder, the hard-core boson model in two-dimensions maps onto a S = 1/2 
XY model and the world-line algorithm for the S = 1/2 XY model can be mapped onto 



the loop algorithm |14] for the 6- vertex model JT5[. When used for world-lines, this method 
replaces the local movement of the world-lines by global changes of loops: the world-lines 
decompose into loops, each of which can be "flipped" with a probability of 1/2. In the 
presence of disorder the flipping probability for a loop is 

l/[l + exp(|Ar (l-2ni)Vi] (9) 

i e loop 
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where Ar is the imaginary-time spacing characteristic of the world-line method. This flip- 
ping probability and the loop construction pl3| , |I"5f define our algorithm. Both the cluster 
and the conventional method suffer very long equilibration and autocorrelation times in 
strong disorder and at low temperature. These facts are the limiting factors of the finite- 
temperature simulations. 

Because the winding number of each world-line is conserved in the conventional method, 
somewhat complicated procedures are necessary for measuring the superfluid density from 
the winding number variance ||. In the new algorithm, however, winding number W is not 
conserved, so we can directly measure the superfluid density from p s = (W 2 )/4J/5 where 
W is D(r) in (||) evaluated at r = /3. 

We performed simulations for two sets of space-time aspect ratios: one with (3/L 2 = 1/4 
and the other with /3/L 1 / 2 = 1/2, corresponding to the two predictions z = 2 and z = 1/2. 
We used an open boundary condition in the temporal direction in the latter case to make 
it possible for the winding numbers to take non-integral values. Allowing this possibility is 
useful when the winding number fluctuations are much smaller than unity, as in the case 
of the second set of simulations. The number of MC steps used was up to 2 x 10 5 for the 
severest case. The number of random samples for the average is 32 to 128. We carefully 
monitored the equilibration times by calculating a Hamming distance similar to the one 
adopted in [|J and found that as the system size or the inverse temperature increases, the 
equilibration time increases rapidly. 

The finite-size scaling form for the finite-temperature superfluid density is 

Ps (6,L } (3)=[3- 1 r(5L 1 /»,(3/L z ), (10) 

with the scaling function T. The critical value V c should be the same as that obtained in 
the zero-temperature calculations. In Fig. |3|, we present results corresponding to the two 
sets of simulations. Based on flip]) , curves should display a common crossing at V = V c if 
the chosen aspect ratio corresponds to the correct z. 

We find that, contrary to the observation in ||, a common crossing point exists at z = 2. 

7 



In addition, the location of the estimated crossing point in (a) is V c = 9.78 ± 0.35, in good 
agreement with the zero-temperature calculation. An apparent crossing occurs in (b) for 
z — 1/2, but the crossing points have expanded to fill a larger region. Such "false" crossings 
can be difficult to diagnose given a limited range of system sizes and temperatures. Here, 
however, we note that the crossing region in (b) is very different from that found in the 
zero-temperature simulation, which should not be the case if z = 1/2. Finally, we remark 
that an analysis of the results in (a) yields v = 0.90 ± 0.13. 

In conclusion, this combination of zero- and finite-temperature algorithms have proven 
very powerful in studies of the superfluid-insulator transition in disordered hard-core bosons. 
The GFMC algorithm we have used produces exact results for the energy, compressibility, 
and superfluid density. As a zero-temperature method, it allows us to perform scaling in the 
spatial variables only, eliminating the need for any a priori assumptions. Our results for the 
compressibility and superfluid density convincingly demonstrate that z — 2. This conclusion 
is strongly supported by the finite-temperature simulations, which are completely consistent 
with the zero-temperature calculations. The combination of zero- and finite-temperature 
simulations allow for a quite sensitive determination of the correct scaling behavior. Further 
details concerning our algorithms and our simulations will be reported elsewhere. 
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FIGURES 

FIG. 1. The superfluid density p s for different lattice sizes as a function of the strength of the 
disorder. The inset shows the corresponding compressibility k in the vicinity of V c . 

FIG. 2. Scaling plots assuming (a) z = 2, (b) z = 1/2, and (c) z = 3. The common crossing 
point is evident for z = 2. For z = 1/2, the crossing points systematically shift to the right as L is 
increased. For z = 3, the shift is to the left. 



FIG. 3. Scaling plots with fixed aspect ratio assuming (a) z = 2 and (b) z = 1/2. 
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